##############################################################################

library(data.table)
library(argparser)
options (warn = -1)

##############################################################################

argp <- arg_parser("Plot the deconstructSigs")
argp <- add_argument(argp, "--input_file" , help="")
argp <- add_argument(argp, "--mutpos_file" , help="")
argp <- add_argument(argp, "--out_file" , help="")
argp <- add_argument(argp, "--mutType" , help="")
argp <- add_argument(argp, "--mut_length" , help="")

argv <- parse_args(argp)

out_file <- argv$out_file
mutType <- argv$mutType
input_file <- argv$input_file
mutpos_file <- argv$mutpos_file
mut_length <- as.numeric(argv$mut_length)


if(1!=1){

  input_file <- "/public/home/xxf2019/20220205_lungSomatic/results/vcf_qc/tmp/GS6_GS1/GS6_GS1.chrX:150523587-150523587.distance.txt"
  mutpos_file <- "/public/home/xxf2019/20220205_lungSomatic/results/vcf_qc/tmp/GS6_GS1/GS6_GS1.chrX:150523587-150523587.mutpos.txt"
  out_file <- "/public/home/xxf2019/20220205_lungSomatic/results/vcf_qc/tmp/GS6_GS1/GS6_GS1.chrX:150523587-150523587.distance.calculate.txt"
  mutType <- "SNV"

  input_file <- "/public/home/xxf2019/20220205_lungSomatic/results/vcf_qc/tmp/GS6_GS1/GS6_GS1.chr22:47703650-47703650.distance.txt"
  mutpos_file <- "/public/home/xxf2019/20220205_lungSomatic/results/vcf_qc/tmp/GS6_GS1/GS6_GS1.chr22:47703650-47703650.mutpos.txt"
  out_file <- "/public/home/xxf2019/20220205_lungSomatic/results/vcf_qc/tmp/GS6_GS1/GS6_GS1.chr22:47703650-47703650.distance.calculate.txt"
  mutType <- "1I"

  input_file <- "/public/home/xxf2019/20220205_lungSomatic/results/vcf_qc/tmp/GS6_GS1/GS6_GS1.chr1:1196071-1196071.distance.txt"
  mutpos_file <- "/public/home/xxf2019/20220205_lungSomatic/results/vcf_qc/tmp/GS6_GS1/GS6_GS1.chr1:1196071-1196071.mutpos.txt"
  out_file <- "/public/home/xxf2019/20220205_lungSomatic/results/vcf_qc/tmp/GS6_GS1/GS6_GS1.chr1:1196071-1196071.distance.calculate.txt"
  mutType <- "24I"

  input_file <- "/public/home/xxf2019/20220205_lungSomatic/results/vcf_qc/tmp/GS6_GS1/GS6_GS1.chr1:4990618-4990618.distance.txt"
  mutpos_file <- "/public/home/xxf2019/20220205_lungSomatic/results/vcf_qc/tmp/GS6_GS1/GS6_GS1.chr1:4990618-4990618.mutpos.txt"
  out_file <- "/public/home/xxf2019/20220205_lungSomatic/results/vcf_qc/tmp/GS6_GS1/GS6_GS1.chr1:4990618-4990618.distance.calculate.txt"
  mutType <- "1D"

  input_file <- "/public/home/xxf2019/20220205_lungSomatic/results/vcf_qc/tmp/GS6_GS1/GS6_GS1.chr1:1196066-1196066.distance.txt"
  mutpos_file <- "/public/home/xxf2019/20220205_lungSomatic/results/vcf_qc/tmp/GS6_GS1/GS6_GS1.chr1:1196066-1196066.mutpos.txt"
  out_file <- "/public/home/xxf2019/20220205_lungSomatic/results/vcf_qc/tmp/GS6_GS1/GS6_GS1.chr1:1196066-1196066.distance.calculate.txt"
  mutType <- "SNV"

  input_file <- "/public/home/xxf2019/20220205_lungSomatic/results/vcf_qc/tmp/GS6_GS1/GS6_GS1.chr1:37299515-37299515.distance.txt"
  mutpos_file <- "/public/home/xxf2019/20220205_lungSomatic/results/vcf_qc/tmp/GS6_GS1/GS6_GS1.chr1:37299515-37299515.mutpos.txt"
  out_file <- "/public/home/xxf2019/20220205_lungSomatic/results/vcf_qc/tmp/GS6_GS1/GS6_GS1.chr1:37299515-37299515.distance.calculate.txt"
  mutType <- "SNV"

  input_file <- "/public/home/xxf2019/20220205_lungSomatic/results/vcf_qc/tmp/GS6_GS1/GS6_GS1.chr1:56954431-56954431.distance.txt"
  mutpos_file <- "/public/home/xxf2019/20220205_lungSomatic/results/vcf_qc/tmp/GS6_GS1/GS6_GS1.chr1:56954431-56954431.mutpos.txt"
  out_file <- "/public/home/xxf2019/20220205_lungSomatic/results/vcf_qc/tmp/GS6_GS1/GS6_GS1.chr1:56954431-56954431.distance.calculate.txt"
  mutType <- "MBS"
  mut_length <- 3

  input_file <- "/public/home/xxf2019/20220205_lungSomatic/results/vcf_qc/tmp/DN1707110057LD01_DN1707110068LD01/DN1707110057LD01_DN1707110068LD01.chr9:92416084-92416084.distance.txt"
  mutpos_file <- "/public/home/xxf2019/20220205_lungSomatic/results/vcf_qc/tmp/DN1707110057LD01_DN1707110068LD01/DN1707110057LD01_DN1707110068LD01.chr9:92416084-92416084.mutpos.txt"
  out_file <- "/public/home/xxf2019/20220205_lungSomatic/results/vcf_qc/tmp/DN1707110057LD01_DN1707110068LD01/DN1707110057LD01_DN1707110068LD01.chr9:92416084-92416084.distance.calculate.txt"
  mutType <- "SNV"  

  input_file <- "/public/home/xxf2019/20220205_lungSomatic/results/vcf_qc/tmp/DN1707110057LD01_DN1707110068LD01/DN1707110057LD01_DN1707110068LD01.chr5:106002624-106002624.distance.txt"
  mutpos_file <- "/public/home/xxf2019/20220205_lungSomatic/results/vcf_qc/tmp/DN1707110057LD01_DN1707110068LD01/DN1707110057LD01_DN1707110068LD01.chr5:106002624-106002624.mutpos.txt"
  out_file <- "/public/home/xxf2019/20220205_lungSomatic/results/vcf_qc/tmp/DN1707110057LD01_DN1707110068LD01/DN1707110057LD01_DN1707110068LD01.chr5:106002624-106002624.distance.calculate.txt"
  mutType <- "1D"
  mut_length <- 1


}

##############################################################################
dat <- data.frame(fread(input_file , header = F))
dat_mutpos <- data.frame(fread(mutpos_file , header = F))

## 只有M认为是比对长度，剪接、插入、缺失都不算在长度里

result_dis <- c()
## 若为单碱基的改变，CIGAR无显示信息
if( mutType=="SNV"){

	for( i in 1:length(dat$V1) ){
		
		tmp <- dat$V1[i]
		mut_pos <- dat_mutpos$V1[i]

		## 计算与参考碱基匹配的长度
		length <- sum(as.numeric(
				sapply(
					strsplit(grep( "M" , unlist(strsplit(tmp , "[0-9]*[H|S]"))  , value = T) , "[M|I|D]") , 
				"[" )
			))

		left_dis <- mut_pos
		right_dis <- length - left_dis

		shortest_dis <- min( left_dis , right_dis )
		result_dis <- c( result_dis , shortest_dis )
	}

}else if( mutType=="MBS"  ){
	## 若为多碱基的替换
	for( i in 1:length(dat$V1) ){
		
		tmp <- dat$V1[i]
		mut_pos <- dat_mutpos$V1[i]

		## 计算与参考碱基匹配的长度
		length <- sum(as.numeric(
				sapply(
					strsplit(grep( "M" , unlist(strsplit(tmp , "[0-9]*[H|S]"))  , value = T) , "[M|I|D]") , 
				"[" )
			))

		left_dis <- mut_pos
		right_dis <- length - left_dis - mut_length + 1

		shortest_dis <- min( left_dis , right_dis )
		result_dis <- c( result_dis , shortest_dis )
	}
}else{
	## 若为插入缺失，CIGAR包含信息
	for( tmp in dat$V1 ){

		## 原先的版本，当一个read上有两个相同的插入和缺失
		if(1!=1){
				## 按照CIGAR分开，分别计算到两边的距离
				left_dis <- unlist(strsplit(tmp , mutType))[1]
				right_dis <- unlist(strsplit(tmp , mutType))[2:length()]

				## 计算与参考碱基匹配的长度
				left_dis <- sum(as.numeric(
						sapply(
							strsplit(grep( "M" , unlist(strsplit(left_dis , "[0-9]*[H|S]"))  , value = T) , "[M|I|D]") , 
						"[" )
					))

				right_dis <- sum(as.numeric(
						sapply(
							strsplit(grep( "M" , unlist(strsplit(right_dis , "[0-9]*[H|S]"))  , value = T) , "[M|I|D]") , 
						"[" )
					))

				shortest_dis <- min( left_dis , right_dis )
				result_dis <- c( result_dis , shortest_dis )
		}
	}

	## 基于突变的位置，考虑碱基的长度
	for( i in 1:length(dat$V1) ){
		
		tmp <- dat$V1[i]
		mut_pos <- dat_mutpos$V1[i]

		## 计算与参考碱基匹配的长度
		length <- sum(as.numeric(
				sapply(
					strsplit(grep( "M" , unlist(strsplit(tmp , "[0-9]*[H|S]"))  , value = T) , "[M|I|D]") , 
				"[" )
			))

		left_dis <- mut_pos
		right_dis <- length - left_dis - mut_length + 1

		shortest_dis <- min( left_dis , right_dis )
		result_dis <- c( result_dis , shortest_dis )
	}
}


## https://www.demo2s.com/java/apache-commons-descriptivestatistics-getelement-int-index.html

DistanceToAlignmentEndMedian <- median(result_dis)
DistanceToAlignmentEndMAD <- median( abs(result_dis - DistanceToAlignmentEndMedian))

result <- data.frame( DistanceToAlignmentEndMedian = DistanceToAlignmentEndMedian , DistanceToAlignmentEndMAD = DistanceToAlignmentEndMAD )

##############################################################################
write.table( result , out_file , row.names = F , sep = "\t" , quote = F )



